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Quadratic-response theory is shown to provide a conceptually simple but accurate approximation 
£ — ■ for the self-consistent one-electron potential of semiconductor nanostructures. Numerical examples 

are presented for GaAs/AlAs and Ino.53Gao.47As/InP (001) superlattices using the local-density 
approximation to density-functional theory and norm-conserving pseudopotentials without spin- 
£N) ■ orbit coupling. When the reference crystal is chosen to be the virtual-crystal average of the two bulk 

constituents, the absolute error in the quadratic-response potential for Fis valence electrons is about 
2 meV for GaAs/AlAs and 5 meV for Ino.53Gao.47As/InP. Low-order multipole expansions of the 
electron density and potential response are shown to be accurate throughout a small neighborhood 
of each reciprocal lattice vector, thus providing a further simplification that is confirmed to be valid 
for slowly varying envelope functions. Although the linear response is about an order of magnitude 
larger than the quadratic response, the quadratic terms are important both quantitatively (if an 
accuracy of better than a few tens of meV is desired) and qualitatively (due to their different 
C/3 ' symmetry and long-range dipole effects). 

PACS numbers: 73.21.-b, 73.61.Ey, 71.15.Ap 
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I. INTRODUCTION 

The potential energy of the electron in the single- 
electron approximation provides the theoretical founda- 
tion for most modeling and interpretation of the elec- 
tronic energy bands of nonmagnetic semiconductors. 
Implementations of this concept range in sophistica- 
tion from simple empirical pseudopotentials^ to first- 
principles quasiparticle self-energies^^ Modern research 
on semiconductor physics deals mainly with complex sys- 
tems (such as nanostructures and novel materials), al- 
though the spintronic properties of traditional semicon- 
ductors are also of strong current interest. It is therefore 
somewhat surprising that, in this day of supposedly ad- 
vanced understanding, no clear physical picture exists of 
something so basic as the one-electron potential energy 
at an ideal GaAs/AlAs or GalnAs/InP heterojunction. 

The key issue is the self-consistent redistribution of 
charge that must occur at any interface between two bulk 
semiconductors. First-principles numerical calculations 
based on density-functional theory and norm-conserving 
pseudopotentials^ are now routinely performed on small 
systems (say, a few tens of atoms) using personal comput- 
ers and free software However, such computations 
are still not feasible for large quantum dots and quantum 
wires. Even if they were, there is no means at present 
for separating this information from a huge numerical 
calculation. That is, no one has yet demonstrated an 
accurate method for breaking down the self-consistent 
potential into simple conceptual units that can be used 
to understand the interface potential and model it us- 
ing less numerically intensive methods (such as empirical 
pseudopotentials or envelope- function theory). 

The purpose of this paper is to demonstrate that the 
quadratic-response formalism developed in Refs. [13 and 
[ill provides a conceptually simple method for construct- 
ing the interface potential that is accurate to within a 



few meV in typical semiconductor heterostructures. In 
this approach, the heterostructure is treated as a pertur- 
bation of a bulk reference crystal, and the self-consistent 
heterostructure potential is built up as a series of succes- 
sively smaller corrections to the bulk potential. 

The leading term is the linear response, which has the 
form of a superposition of elementary quasi-atomic build- 
ing blocks. This form is widely used in empirical pseu- 
dopotential models (see the literature review in Sec. HI)) , 
but it is shown here to be accurate to only a few tens 
of meV. The quadratic response is then added in the 
form of diatomic building blocks, which are qualitatively 
different because they carry a dipole moment even in cu- 
bic semiconductors. Numerical implementations of the 
theory are presented here for the paradigmatic common- 
atom and no-common-atom systems GaAs/AlAs and 
Ino.53Gao.47As/InP. The quadratic-response approxima- 
tion for Ti5 valence electrons is shown to be accurate 
to within about 2 meV for GaAs/AlAs and 5 meV for 
Ino.53Gao.47As/InP. 

In typical heterostructures, the low-energy electron 
and hole states can be represented by slowly varying 
envelope functions. The properties of such states are 
determined by the pseudopotential not just at the bulk 
reciprocal-lattice vectors G (as would be the case in a 
bulk semiconductor), but throughout a small k-space 
neighborhood of each G. Within such a neighborhood, 
an analytic function of k can be approximated using a 
power series (which is equivalent to introducing a mul- 
tipole expansion of a localized function in coordinate 
space). This paper also examines the accuracy of such 
power series and demonstrates that the truncated expan- 
sions introduced in Refs. 03 and[ll] are sufficiently accu- 
rate for slowly varying envelope functions. The appli- 
cation of these multipole expansions in actual envelope- 
function calculations is presented in another paper 

The paper begins in Sec. [Til with a review of previ- 
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ous models of the interface potential. Section IIIII de- 
scribes the choices made in defining the model system 
used here for numerical calculations. The basic features 
of the self-consistent density and potential in GaAs/AlAs 
and Ino.53Gao.47As/InP superlattices are examined in 
Sec. IIV( here the valence-band offsets calculated for the 
model systems are found to be in reasonably good agree- 
ment with experiment. Section fVl outlines the methods 
used to calculate the linear and quadratic response to the 
hetcrostructure perturbation and discusses the physical 
significance of the results. Localization of the density and 
potential response and the use of multipole expansions 
at small wave vectors are considered in Sec. I VI I Section 
IVIII shows that the cumulative errors arising from the 
quadratic-response and truncated-multipole approxima- 
tions are limited to a few meV for the systems studied. 
The conclusions of the paper are reviewed in Sec. IVIIII 
Finally, Sec. IIXI discusses several limitations of the calcu- 
lations presented here, along with possible extensions of 
this work that could be used to overcome some of these 
limitations. 



II. PREVIOUS DESCRIPTIONS OF THE 
INTERFACE POTENTIAL 

Two simple non-self-consistent models of the hetero- 
structure pseudopotential are widely used. The simplest 
is the planar interface model, in which the potential is 
assumed to make an abrupt transition from one periodic 
bulk potential to another at some predefined interface 
planej 13 i 14 ' 15 i 16 ' 17 i 18 ' 19 ' 20 This requires as input only the 
bulk pseudopotential form factors?^ The assumption of 
a planar interface provides, in essence, a simple recipe for 
interpolating the pseudopotential in k space between the 
bulk reciprocal-lattice vectors G. A discontinuous po- 
tential is, however, clearly far from self-consistent (since 
by Poisson's equation it would require an infinite elec- 
tron density), and a planar interface fails to account for 
the basic physics of the bonds between atoms at a no- 
common-atom interfaced 

An arguably more realistic potential (although not 
all would agreed) is provided by the atomic superposi- 
tion model, which is defined as a linear combination of 
screened atomic pseudopotentialSi 22 ' 23 i 24 ' 25 i 26 ' 27 i 28 Here 
the atomic potentials are assumed to be known as contin- 
uous functions of kj 2 ^ hence, no interpolation is required. 
In its simplest form this model may not accurately re- 
produce differences between bulk semiconductors, 26 since 
(for example) the screened As pseudopotential obtained 
by fitting the properties of bulk GaAs will generally differ 
from that obtained by fitting AlAs. 

However, Mader and Zunger have overcome this dif- 
ficulty by using environment-dependent pseudopoten- 
tials^ 6 - in which the replacement of one Ga atom in GaAs 
by Al is accompanied by small changes in each of the four 
neighboring As pseudopotentials. The atomic pseudopo- 
tential and its associated nearest-neighbor perturbations 



may be grouped conceptually into a single quasi-atomic 
unit, which has the site symmetry Td of an Al impurity 
in GaAs rather than the spherical symmetry of an iso- 
lated Al atom. Thus, in this more general quasi-atomic 
superposition model^ 6 - the heterostructure potential is 
constructed as a linear combination of quasi-atomic po- 
tentials, each of which has the site symmetry of a substi- 
tutional impurity. (The planar interface model may also 
be interpreted as a superposition of elementary units, 30 
but in this case the units do not have the atomic site 
symmetry.) 

The superposition potential is also manifestly non-self- 
consistent, if only because the exchange-correlation po- 
tential is nonlinear. The assumption of linearity has im- 
portant physical consequences. As shown in Eq. (10) 
of Ref. [3l|, one of these is the absence of direct Ti- 
X\ intervalley coupling in GaAs/AlAs (001) heterostruc- 
tures. The same conclusion was reached in Ref. using 
a superposition of quasi-atomic potentials that were con- 
structed by a different method. These conclusions were 
criticized by Takhtamirov and VolkoA*^ 3 - on the grounds 
that the true interface potential does not have the as- 
sumed form, and that there is consequently no basis for 
assuming T\-Xi coupling to be any smaller than Ti 
A3 coupling (the latter of which is nonvanishing in the 
quasi-atomic superposition mode l 31 i 32 ). This criticism is, 
of course, entirely correct if the self-consistent interface 
potential is viewed as completely unknown. 

Indeed, many derivations of envelope-function mod- 
els are prefaced by a warning that the potential within 
a few lattice constants of the interface is too compli- 
cated for direct analysiS ) 34 i 35 ' 36 i 37 ' 38 i 39 even in the case 
of an ideal interface (i.e., no interface roughness, de- 
fects, or interdiffusion) . In early works it was some- 
times concluded that the details of the interface potential 
have no significant influence on slowly varying envelope 
functions. Today, however, it is widely accepted that 
these details are of critical importance in determining 
the strength of small symmetry-breaking interface terms 
in the Hamiltonian, such as the above-mentioned in- 
tervalley mixing ) 31 i 32 i 33 valence-band mixing ) 32 i 40 ' 41 and 
conduction- and valence-band Rashba coupling 4 ^ (see the 
following paper— for further study of some of these top- 
ics). But the assumption that the self-consistent poten- 
tial cannot be understood in simple terms remains un- 
challenged by most envelope-function theorists. 

This problem was tackled in a landmark se- 
ries of papers by Baroni et al, who showed that 
linear-response theory is very successful in predict- 
ing valence-band offsets in semiconductor heterostruc- 
turesM^d^dL^dl^l^ i n this approach, the ionic 
pseudopotential of the heterostructure is treated as a per- 
turbation of a bulk reference crystal. The linear response 
to this perturbation was shown to give accurate predic- 
tions for the valence-band offset (within about 0.01 eV 
of the value obtained from a direct calculation, and in 
good agreement with experiment) if the reference crystal 
was chosen to be the virtual-crystal average of the two 
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bulk materials, since the quadratic response is then the 
same on both sides of the heterojunction and gives no 
contribution to the offset. However, if one is interested 
in the actual position-dependent potential, including the 
details of its behavior near the interface, the quadratic 
response is not generally negligible. 

Nonlinear screened empirical pseudopotentials were 
introduced by Magri and Zunger,— who extended the 
environment-dependent concept of Mader and Zunger 26 
to the no-common-atom InAs/GaSb system. Nonlinear- 
ity arises in this case because the nearest-neighbor per- 
turbations depend on which atoms occupy the neighbor- 
ing sites. The pseudopotentials can therefore no longer 
be interpreted as purely quasi-atomic units; a complete 
description requires the inclusion of nearest-neighbor di- 
atomic building blocks as welli^ This represents an im- 
portant conceptual advance in the modeling of hetero- 
structure pseudopotentials. 

The present work develops this concept more thor- 
oughly by extending the first-principles linear-response 
theory of Baroni et a/.— to include the diatomic quadratic 
response for arbitrary pairs of atoms. The results show 
that the environment dep endence of the empirical pseu- 
dopotentials in Refs. l26l and lo3l is qualitatively incom- 
plete, since (for example) the leading contribution to T\- 
X\ mixing in GaAs/AlAs (001) superlattices arises from 
second-nearest-neighbor Ga-Al pairs. (However, this 
contribution is about two orders of magnitude smaller 
than the linear response, which supports the conclusion 
that Y\—X\ mixing is much weaker than Y\-X-s mix- 
ing.) Also, the diatomic nearest-neighbor response for 
no-common-atom systems calculated here has a dipole 
moment that generates long-range electric fields extend- 
ing over the entire superlattice. Such fields do not oc- 
cur in the empirical pseudopotential interpolation scheme 
of Magri and Zunger ,~ which produces only short-range 
potentials. This suggests that a more complete descrip- 
tion of quadratic diatomic response may be of benefit in 
achieving more realistic empirical pseudopotentials. 



III. NUMERICAL CONSIDERATIONS 
A. Definition of model system 

All of the numerical results in this paper are 
derived from plane-wave pseudopotential total-energy 
calculations 5,6 performed using the ABINIT software.^^ 
This package provides a variety of options, but the par- 
ticular calculations reported here used the local den- 
sity approximation (LDA) to density-functional theory 
and the nonlocal norm-conserving pseudopotentials of 
Hartwigsen, Goedecker, and Hutter— with no spin-orbit 
coupling. 

Choosing a particular physical structure and a par- 
ticular set of technical ingredients (such as plane-wave 
kinetic-energy cutoffs and k-point sampling) defines a 
model system whose properties can be calculated self- 
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Kinetic energy (hartrees) 



FIG. 1: (Color online) Kinetic energy divisor (with continu- 
ous slope) and nonlocal potential energy multiplier (with con- 
tinuous curvature) used to smooth out discontinuities in the 
momentum-space electron density and pseudopotential. The 
multiplier used for the local pseudopotential has the same 
shape, but the energy scale is four times that shown. 



consistently to an accuracy approaching that of machine 
precision. These "exact" model calculations are used 
here and in the following paper— as a benchmark for 
comparison with various approximations used in the con- 
struction of the first-principles envelope-function theory 
of Ref. [Til. The main objective of this paper is to pin 
down as closely as possible how much error arises from 
the quadratic-response approximation and the multipole 
expansions of the linear and quadratic response. The ap- 
proximate potentials derived here are also used directly 
in a subsequent paper— that examines the accuracy of 
first-principles envelope-function theory in superlattices. 

With these goals in mind, relatively low kinetic-energy 
cutoffs (see below) and numbers of special k points were 
chosen in order to make the "exact" calculations feasible 
for the large superlattices in which envelope- function the- 
ory is valid. In particular, two independent k points were 
used for GaAs/AlAs superlattices with Did symmetry, 
while four points were used for Ino.53Gao.47As/InP super- 
lattices with Civ symmetry. The total energy is not com- 
pletely converged at these values, but the valence-band 
offsets calculated here are nevertheless in good agreement 
with experiment and with previous first-principles calcu- 
lations reported in the literature. 

In order to calculate multipole moments of the linear 
and quadratic response, it is necessary that the electron 
density and short-range potential response be well local- 
ized within the given large supercells. However, sharp 
cutoffs in the plane-wave basis produce spurious long- 
range Gibbs oscillations that are very small but nonethe- 
less large enough to completely swamp the physical mul- 
tipole values. To suppress these oscillations, the smooth- 
ing functions shown in Fig. [TJ were applied to the kinetic 
and potential energy. The chosen kinetic-energy cutoff 
(i.e., the value at which the divisor in Fig.[T]goes to zero) 
includes 283 plane waves at the T point in bulk material 
(113 of which are not altered by smoothing) and more 
than 14000 plane waves in a 96-atom superlattice. 
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B. Conventions denning the absolute energy 

There are two conventions in the momentum-space 
total-energy formalism that determine the absolute value 
of the local potential. First, the mean value of the 
Hartree potential is set to zero by definition i 56 i 57 i 58 

Wq) = (l-5 q o)^n(q) (1) 

(where n is the valence electron density), because in a 
neutral system the net contribution to the total energy 
from the Coulomb interaction terms (for the electron-ion 
system) at wave vector q = is identically zero. 

A similar convention 5 -^ is often used for the mean value 
of the short-range part of the local pseudopotential. This 
is given by V psp = J2j Vj > where^&S&iSi 

Vj = hJ (w r )+f) d>r ' ( 2 ) 

in which f2 is the unit cell volume and Vj and Zj are 
the local pseudopotential and valence charge of ion j. 
The contribution from V^ sp to the total energy is just 
the constant -E^orc = ZQV psp , where Z = Y) j Zj. Hence, 
it is permissible to set V pBp to zero in the Kohn-Sham 
equations^ 

However, in order to calculate the nonlinear response 
correctly, V psp must be included in the local pseudopoten- 
tial. In the present work, V psp was added to both the lo- 
cal pseudopotential and the energy eigenvalues after the 
conclusion of each self-consistent calculation. Therefore, 
the absolute energy scale for all figures in this paper (and 
the nexlji^) is defined solely by the convention of setting 
the mean Hartree potential to zero. 



IV. VALENCE-BAND OFFSETS IN 
SUPERLATTICES 

Energy band structures for the model system described 
in Sec. Mil are presented in Figs. [2] and [3] for the bulk ma- 
terials GaAs, AlAs, Ino.53Gao.47As, and InP. The con- 
duction bands in these figures are clearly not realistic, 
as the minimum at L lies slightly below the T minimum 
in GaAs, and is only slightly above it in Ino.53Gao.47As 
and InP. The calculated ri c -Ti5 V energy gaps are 1.850 
eV for GaAs and 1.492 eV for Ino.53Gao.47As, which are 
13% and 61% larger than the experimental values 5 -^ (with 
spin-orbit splitting removed) of 1.633 eV and 0.925 eV. 
Nevertheless, the model does capture many of the gross 
qualitative trends in experimentally determined band 
structures^ such as the predominance of the X valley 
in AlAs relative to GaAs. 

It is important to note that the difference between the 
valence-band maxima in Figs. [2] and [3] is not equal to the 
valence-band offset^ 5 - To obtain the valence-band offset, 
one must add to this difference the shift in the macro- 
scopic Hartree potential of the two materials as given by 
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FIG. 2: (Color online) Energy band structure of bulk GaAs 
and AlAs. 
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FIG. 3: (Color online) Energy band structure of bulk 
Ino.53Gao.47As and InP. 

a supercell calculation. 45 The latter shift is shown in Figs. 
|4] and [5] for a GaAs/ AlAs superlattice and in Figs. [6] and 
[7] for an Ino.53Gao.47As/InP superlattice. These figures 
show planar averages and macroscopic average o 60 ' 61 i 62 of 
the electron density and Hartree potential, with distance 
given in units of the monolayer spacing d — ^a, where a 
is the cubic lattice constant. 

To calculate the macroscopic average, the microscopic 
functions were first averaged over the volume of a bulk 
unit cell, as in Ref. [6(1 Then, to ensure that the out- 
put is smooth and finite even when the input is singular 
(e.g., in the case of a truncated multipole expansion), an 
additional averaging was performed with respect to the 
normalized Gaussian function 61 

/(*)= ±exp[-7r(z/<7) 2 ], (3a) 
a 

which is equivalent to multiplying in momentum space 
by 

g(k) = exp[-7r(/ccr/27r) 2 ]. (3b) 
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FIG. 4: (Color online) Density and potential of a (001) FIG - 6: ( Color online ) Density and potential of a (001) 

(GaAs)24(AlAs) 2 4 superlattice: (a) Planar average of elec- (Ino.53Gao.47As)24(InP) 2 4 superlattice: (a) Planar average of 

tron density, (b) Planar average of Hartree potential, (c) electron density, (b) Planar average of Hartree potential, (c) 

Macroscopic average of density and potential. Macroscopic average of density and potential. 




FIG. 5: (Color online) Expanded view of the planar average 
density and potential from Fig. [4] 



In Figs. |4jc) and [6^c) the distance a was chosen to be 
a = d. 

The macroscopic density for GaAs/AlAs in Fig. [4jc) 
has the familiar interface dipole shape,— ^ leading to a 
steplike shift in the Hartree potential of 0.412 eV. The 
shift predicted by linear response theory is very close to 



this value, at 0.416 eV (see below for a more detailed 
study of the error in this approximation) . These results 
agree well with the shift of 0.41 eV calculated in Ref. liH . 
Since the difference between the valence-band maxima 
in Fig. [2] is 18 meV, the valence-band offset for a (001) 
GaAs/AlAs heterojunction is 0.430 eV for the model sys- 
tem, as compared to the value of 0.45 eV reported in Ref. 
|45| . If this figure is corrected for spin-orbit coupling and 
quasiparticle effects using the values in Ref. |4f|, the net 
offset is roughly 0.56 eV, in good agreement with the 
experimental value 59 of 0.53 eV. 

The macroscopic density of the no-common-atom 
Ino.53Gao.47As/InP superlattice in Fig. [6jc) does not 
have a simple dipole shape because (to leading order) it 
is a superposition of two offset dipoles, one for the cation 
and one for the anion i 45 ' 46 As discussed below, there is 
also a slight asymmetry of the two interfaces, leading to 
a macroscopic electric field that is barely visible in Fig. 
Etc). The linear Hartree shifts [from Eq. (O with I = 0] 
for cations and anions are +0.366 eV and —0.697 eV, as 
compared to the values of +0.34 eV and —0.58 eV re- 
ported in Refs.[45l andlirl The net linear Hartree shift of 
—0.331 eV agrees well with the mean shift of —0.334 eV 
in Fig.[5Jc). Combining this with the 0.574 eV difference 
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FIG. 7: (Color online) Expanded view of the planar average 
density and potential from Fig. Here X represents the 
virtual atom Ino.53Gao.47. 
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FIG. 8: (Color online) Linear and quadratic density response 
to a monatomic Ga perturbation of Alo.5Gao.5As. 



in the valence-band maxima of Fig. [3J the net valence- 
band offset for the model Ino.53Gao.47As/InP system is 
0.240 eV. After adjusting for experimental spin-orbit 
splitting^ the calculated offset of 0.313 eV is close to 
the experimental value 59 of 0.345 eV. 



V. LINEAR AND QUADRATIC RESPONSE 

Having established that the model system adopted 
here provides at least a rough approximation of phys- 
ical reality, the next step is to calculate the linear 
and quadratic response to virtual-crystal perturbations 
of a bulk reference crystal. The reference crystal is 
chosen here to be the virtual-crystal average of the 
bulk constituents (i.e., Alo.5Gao.5As for GaAs/AlAs and 
Ino.765Gao.235Aso.5Po. 5 for lno.53Gao.47As/InP). The lat- 
tice constant for all calculations was fixed at the value 
obtained by minimizing the total energy of the reference 
crystal (i.e., a = 10.5 bohr for Alo.5Gao.5As and a = 10.9 
bohr for In .765Ga .235Aso.5Po.5)- 

The perturbation of the heterostructure relative to the 
reference crystal is defined by the change in pseudopo- 
tential 



psp 



'R^ion 



(x - R Q ), 



(4) 



R 



which is written as a local potential for simplicity (see 
Ref. [ll| for the nonlocal case). Here ^^(x) is the ionic 
pseudopotential of atom a, R Q is the position of atom a 
in unit cell R, and 9^ is the change in fractional weight 
of atom a in cell R of the heterostructure relative to 
the reference crystal. Since the total change in fractional 
weight at each site must add to zero, this can be rewritten 
asii 



Ay psp (x) = ^'5> R A< n (x - R Q ), 



(5) 



where the sum covers only independent values of a (e.g., 
either Ga or Al — but not both — in GaAs/AlAs), and 
At!;" n is the change in v" on relative to the reference crys- 
tal. 

If n(x) is the exact density of the heterostructure, the 
linear and quadratic response to virtual perturbations 9^ 
are defined as the derivatives 



A«r(x) 



9n(x) 
d9^ ' 



An^ R ,(x) = i d 2n(x) ,. 
RR 1 ' 2 cX9 R <9<, 



(6a) 
(6b) 



The total density may then be reconstructed from the 
power series 

n(x) = n (0) (x) + n {1) (x) + n (2) (x) + • • • , (7a) 
nW(x)=^«(x), (7b) 



Ct.R 



,(2) 



(7c) 



,R q',R' 



R 



where n(°)(x) is the density of the reference crystal. In 
the present work, this power series was truncated at the 
second order. The derivatives © were calculated by the 
direct supercell method^ in which (for example) the 
Alo.5Gao.5As reference crystal was perturbed by replac- 
ing one or two Alo.sGao.s atoms with either Alo.55Gao.45 
or Alo.45Gao.55. See the Appendix for further details. 

The linear and quadratic density response to a 
monatomic perturbation in Alo.5Gao.5As are shown in 
Fig. [5J The functions plotted here are planar averages 
of An^(x) and An^(x) for the case a ~ Ga. Since 
the supercell used here is the same as that of the (001) 
superlattice in Figs. [4] and [5j the perturbation actually 
consists of a plane of Ga atoms (or rather an infinites- 
imal perturbation toward Ga). This perturbation has 
D2d symmetry, so the planar average shown in Fig. [5] is 
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FIG. 9: Quadratic density response to a diatomic Ga-Ga per- 
turbation of Alo.5Gao.5As. The perturbations are applied at 
z = —0.5 and z = +0.5. 



symmetric. Note that the monatomic quadratic response 
is about an order of magnitude smaller than the linear 
response. 

The quadratic response to a diatomic perturbation in 
Alo.5Gao.5As is shown in Fig. [§J This perturbation also 
has D2d symmetry, so its qualitative features are similar 
to those of the monatomic response (except that it is 
centered on an anion plane rather than a cation plane). 
The perturbation planes in Fig. [5] are separated by one 
monolayer (i.e., \a)- The calculations performed in this 
paper include diatomic perturbations out to a separation 
of four monolayers. Perturbations beyond this distance 
were neglected because they yield contributions of less 
than 0.1 meV. 

The second-nearest-neighbor diatomic response shown 
in Figgis the leading contribution to T\—X\ coupling in 
GaAs/AlAs (001) superlattices. A comparison with Fig. 
[5] shows that this term is nearly two orders of magnitude 
smaller than the linear response, which is the leading 
contribution to Y\- X3 coupling. Hence, these results in- 
dicate that Yx—X\ coupling is indeed much smaller than 
T\-X$ coupling, as suggested in Refs.[3l| and l32l. 

Qualitatively different results are obtained for 
the perturbation of two different atoms of the 
Ino.765Gao.235Aso.5Po. 5 reference crystal in Fig. [TUJ This 
diatomic perturbation has symmetry C^, so its planar 
average is not symmetric. It therefore has a nonvanishing 
dipole moment, which is readily apparent from the fig- 
ure. The dipole generates a rapid change in the Hartree 
potential, as shown in Fig. fTOTb). In order to satisfy the 
periodic boundary conditions, this steplike change in po- 
tential must be accompanied by a macroscopic electric 
field extending over the entire supercel L 63 ' 64 This field 
polarizes the reference crystal, producing small oscilla- 
tions in the density and potential^ that are visible away 
from z — 0. Figure [10] shows calculations performed on 
two supercells, one with 48 monolayers (the same as the 
superlattice in Figs. [H] and [7]) and one with 12 monolayers. 
It is evident that the polarization amplitude is inversely 
proportional to the period of the supercelL 

It is these dipole moments in the diatomic quadratic re- 
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FIG. 10: (Color online) Quadratic response to a diatomic 
As-Ga perturbation of Ino.765Gao.235Aso.sPo.5: (a) electron 
density, (b) Hartree potential. The perturbations are applied 
at z = —0.5 (As) and 2 = +0.5 (Ga). Results are given for 
supercells with periods N = 48 and N = 12. 



sponse that are responsible for the slight interface asym- 
metry and macroscopic electric field referred to in the 
discussion of Fi g, p Such effects have been studied pre- 
viously in Refs. I6a7 numerically) andllol (analytically). 



VI. MULTIPOLE MOMENTS AND 
LOCALIZATION 

A. Definition of multipole moments 

As discussed in the Introduction, slowly varying enve- 
lope functions probe the density and potential only in a 
small neighborhood of each reciprocal-lattice vector G. 
It is therefore convenient to introduce long-wavelength 
approximations for the density and potential response 
in the form of power-series expansions of n(k + G) and 
v(k + G) with respect to k. The expansion coefficients 
can be found numerically by several methods. One is 
a simple polynomial fitting of the discrete fast Fourier 
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transform (FFT). Another is to use analytical manipula- 
tions of the FFT to obtain expansion coefficients in the 
form of multipole moments. 

In a superlattice, the allowed values of k + G satisfy 
k = kz, for small |k|, where z is the direction normal to 
the interface plane. Hence 

»(k + G) = J- ™(x)e- l ( k+G >- x (8a) 



OO 

E( 

1=0 



-ik) l m(G), 



(8b) 



where x is a coordinate on the FFT grid, N-p is the num- 
ber of such points, and the expansion coefficients are 



n ' (G) = ^^ A(x)e 



-iG x 



i^[0(k-G)]*n(k), 

k 



(9a) 
(9b) 



where 0(k) is the FFT of z l . The coefficient ni(G) is, 
to within a factor of Z!, just the multipole momen t 45 ! 61 of 
order 2 l of the function n(x)e _lG x . 



B. Order of terms included 

For small k, the power series (l8b|) can be approximated 
accurately by including only a few terms with small I. 
This paper adopts the approximation scheme defined in 
Refs. [l0| and I 111 in which the power series for the linear 
and quadratic potential response i/ 1 ' (k+G) and tA 2 ) (k+ 
G) are terminated at the upper limits 

= 2, = 0. (10a) 

These power series can be used for the pseudopotential 
and exchange-correlation potential response, which are 
analytic functions of k. However, the Hartree potential 
is nonanalytic, so the expansion ([8b| must be applied 
to the density n(k + G) instead. To obtain an accuracy 
equivalent to (|10a|) . the upper limits for the density power 
series are extended when G = 0: 

j£' 2) (G)=#' 3 > +2S GQ , (10b) 

because of the k~ 2 factor in the Hartree potential ([1]). 

When G^O, the Hartree potential is an analytic func- 
tion of k (for small k), and the power series expansion of 
|k + G|~ 2 generates the following Hartree corrections to 
the analytic potential coefficients: 

-in 
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When G — 0, analytic potential contributions are also 
obtained from the terms 



A«,(0) = -47rn i+2 (0) (I > 0), 



(12) 



where < I < 2 for the linear response and I = for the 
quadratic response. In a superlattice, the only nonana- 
lytic terms arise from the density monopole, dipole, and 
quadrupole coefficients at G — 0. The quadrupole term is 
given by 47rn2(0)<5ko, which merely contributes an overall 
constant shift in potential. The monopole term is zero 
for the isovalent perturbations considered here, while the 
linear dipole term vanishes for perturbations with T<z or 
D2d symmetry. The quadratic dipole term is, however, 
nonzero for diatomic perturbations with Civ symmetry, 
as shown above in Fig. [17)1 



C. Localization 

Since the FFT grid is finite, the multipole moment 
n;(G) always has a well defined value. Nevertheless, if 
the response n(x) is not well localized, the value of ni(G), 
although finite, does not converge to a meaningful value 
in the limit of large supercells. It is therefore necessary 
to check in each case whether the localization is sufficient 
to use Eq. ^ for the desired values of /. 

The localization of the response is not merely a nu- 
merical detail; it is fundamental to the physics of the in- 
sulating statei 67 ' 68 Phenomenological envelope-function 
models typically assume that the heterostructure Hamil- 
tonian can be expressed as a power series in k, but 
the validity of such an assumption depends crucially on 
the localization of the response in an insulatorj 10 ' 11 If 
the response is not localized, the entire foundation for 
the power-series expansion of the Hamiltonian breaks 
down. Establishing that the response is localized is there- 
fore a key element in the justification of the standard 
phenomenological approach to envelope-function theory. 
Numerical error inevitably masks the underlying local- 
ization at some point, but it is important to ensure that 
this error is insignificant for all practical purposes. 

This is done in Fig. [TT] for the linear density response 
of Alo.5Gao.5As in the cases I = 0,2, and 4. The response 
is quasi-exponentially localized over nearly 10 orders of 
magnitude, to a distance of about 17 monolayers from 
the perturbation. At this point it flattens out due to the 
nonanalyticity of the potential-energy smoothing func- 
tion in Fig.m If the potential energy is not smoothed, the 
quasi-exponential localization extends for another order 
of magnitude, at which point numerical error takes over. 
However, the localization obtained here is clearly suffi- 
cient to calculate accurate values of the linear quadrupole 
and hexadecapole moments. On the other hand, if the 
plane-wave cutoff in the kinetic energy is not smoothed, 
z i n^ 1 \z) is not adequately localized, and the hexade- 
capole moment (I — 4) has a nonsensical value for the 
supercell considered here. 
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FIG. 11: (Color online) Log-scale plot of linear density re- 
sponse in Fig. [8] 
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FIG. 12: (Color online) Log-scale plot of linear and quadratic 
pseudopotential plus exchange-correlation potential for a 
monatomic Ga perturbation of Alo.5Gao.5As. 



The linear and quadratic terms in those parts of the 
local potential response that are expected to be local- 
ized (i.e., the ionic pseudopotential and the exchange- 
correlation potential) are shown in Fig. [T2j Here the 
features seen in the linear response for \z\ > 10 are the 
remnants of Gibbs oscillations in the local ionic pseudo- 
potential that remain even after smoothing of the plane- 
wave cutoff. (These Gibbs oscillations are not seen in the 
quadratic response because the ionic pseudopotential is 
purely linear.) The localization here is not quite as good 
as for the density, but it is still sufficient to calculate the 
desired multipole moments. 

The log-scale plot in Fig. [T2] of the quadratic density 
response for C%v diatomic perturbations has additional 
features of interest. The response is not well localized, 
but this is a physical effect arising from the periodic os- 
cillations shown in Fig. 1101 The periodic part can be 
separated from the rest by evaluating it in the unit cell 
farthest from the origin. Subtracting the periodic part 
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FIG. 13: (Color online) Log-scale plot of the diatomic density 
response in Fig. llOf a). Both the total density and its localized 
part (excluding the periodic part) are shown. 
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FIG. 14: (Color online) Linear and quadratic density response 
to a monatomic Ga perturbation of Alo.sGao.s As. The solid 
lines are Fourier transforms of the functions in Fig. [8] The 
dotted and dashed lines are quadratic and quartic polynomials 
defined by the quadrupole and hexadecapole moments. 



from the total leaves a remainder that should be local- 
ized. This is confirmed in Fig. 1131 which shows both the 
total quadratic response and the localized remainder. 

The localized part of the diatomic response can be used 
to obtain multipole moments as above. The periodic part 
is even simpler, because it merely modifies the values 
of coefficients in the bulk Hamiltonian (as discussed in 
greater detail below). 



D. Calculated multipole moments 

The results obtained from the calculated density mul- 
tipole moments n/(G = 0) are shown in Figs. [TH and 
[13 for Alo.5Gao.5As and ln .765Ga .235Aso.5Po.5, respec- 
tively. The truncated power series expansions for the 
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FIG. 15: (Color online) Fourier transform of the diatomic 
density response in Fig. UOf a). Solid lines are for a super- 
cell with N = 48 monolayers, while circles are for N = 12. 
The FFT is scaled by N so that the results are compara- 
ble. Dashed lines are quadratic (real) and linear (imaginary) 
polynomials corresponding to the calculated quadrupole and 
dipole moments. 



k-space density response are clearly valid only for small 
wave vectors. For the C-i v perturbation in Fig. 1151 the 
reduced symmetry generates a nonvanishing imaginary 
part in the FFT. This part has a nonzero slope at the 
origin, leading to the dipole moment observed previously 
in Fig. [TO] 



There are also noticeable spikes in the FFT at nonzero 
reciprocal-lattice vectors. These correspond to the pe- 
riodic part of the response discussed above (Sec. IVI C"j) . 
whose physical origin is the polarization of the reference 
crystal by the macroscopic electric field^ in Fig. |TD] The 
amplitude of the spikes is obtained most accurately by 
the method described above, in which the periodic part 
of the quadratic response is evaluated in coordinate space 
in the unit cell farthest from the origin and then Fourier 
transformed. The resulting Fourier series coefficients oc- 
cur only at the bulk reciprocal-lattice vectors G, which 
means that their contribution to the envelope-function 
Hamiltoniar>i2, has the same formal structure as the bulk 
kp Hamiltonian of the reference crystal. However, these 
bulk polarization terms have the Ci v symmetry of the 
dipole, not the Td symmetry of the reference crystal. As 
shown in Ref. [HI, the net effect of the Ci v dipoles and 
their associated bulk polarization is to contribute about 
one third of the total splitting of the quasidegenerate X 
and Y valence states (which would be degenerate under 
D2d symmetry) in Ino.53Gao.47As/InP superlattices. 
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FIG. 16: (Color online) Errors in the electron density and 
Hartree potential of a GaAs/AlAs superlattice under various 
approximations: (a) Errors in the reference crystal density 
and the densities obtained in the linear and quadratic ap- 
proximations, (b) Expanded view of the linear and quadratic 
errors in (a), (c) Errors in the linear and quadratic approxi- 
mations to the Hartree potential. 



VII. ERROR IN THE LINEAR AND 
QUADRATIC APPROXIMATIONS 

A. Electron density and Hartree potential 

To test the validity of truncating the power series in 
the nonlinear response ([7]), the effects of various trun- 
cations were examined by comparing them with the ex- 
act superlattice density and potential given previously in 
Figs. [H [5J [H and [7] The approximations are generally 
good enough that they are difficult to distinguish visually 
from the exact results. Therefore, only the error in the 
approximations is plotted here. 

The results for the GaAs/AlAs superlattice are shown 
in Fig. [TH1 Part (a) shows the error when the super- 
lattice density is approximated by the periodic reference 
crystal density, comparing it with the error in the linear 
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FIG. 17: (Color online) Errors in the electron density and 
Hartree potential of an Ino.53Gao.47As/InP superlattice under 
various approximations: (a) Errors in the linear and quadratic 
densities, (b) Errors in the linear and quadratic Hartree po- 
tentials. 
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FIG. 18: (Color online) Sum of the linear and quadratic per- 
turbations for (a) the electron density and (b) the Hartree 
potential constructed from quadrupole and hexadecapole mo- 
ments in a GaAs/AlAs superlattice. The perturbations are 
finite only because a discrete coordinate grid is used. The 
hexadecapole terms are dominant here, although this is no 
longer true after macroscopic averaging. 



and quadratic approximations. The latter are both small 
on the scale of the reference crystal error (which is itself 
small on the scale of the reference density), so the linear 
and quadratic errors are plotted on an expanded scale in 
part (b). The corresponding errors in the Hartree po- 
tential are shown in part (c). It can be seen that the 
quadratic approximation is quite accurate, with an error 
of only about ±2 meV. 

The errors in the linear and quadratic approximations 
for the Ino.53Gao.47As/InP superlattice are shown in Fig. 
[T71 Here the error in the quadratic Hartree potential 
has an oscillation with an amplitude of about 3 meV 
superimposed on a shift of about ±1.5 meV. 



B. Truncation of multipole expansions 

Using multipole expansions for the density and po- 
tential introduces an additional source of error. To il- 
lustrate this, Fig. [TH] shows the electron density and 
Hartree potential in coordinate space calculated directly 
from the FFT of the truncated multipole expansion (with 
< I < 4 for the linear density and < I < 2 for the 
quadratic density). Note that the scale of this figure is 
much larger than any of the others. In fact, in a contin- 
uous coordinate space, these quantities would be deriva- 
tives of Dirac 5 functions. Clearly these multipole ex- 
pansions yield a very poor approximation of the original 



density and potential. 

However, in an envelope-function model, accuracy is 
only required at small wave vectors, so it is more appro- 
priate to compare the macroscopic average of the mul- 
tipole expansions with the macroscopic average of the 
exact results. This is done in Fig. [19l which shows the 
errors in the macroscopic density and Hartree potential 
for a GaAs/AlAs superlattice. Here the macroscopic av- 
erages were calculated using a = 'id in Eq. ([3]), which 
limits momentum transfers to roug hly i of the bulk V- 
X distanced Such a limit is appropriate for the slowly 
varying envelope functions considered in the following 
paper.— The macroscopic errors are virtually identical 
in bulk, but the truncated multipole approximation in- 
troduces some additional error near the interfaces. This 
interface error is negligible for the case a = 3d shown 
here, although it becomes arbitrarily large in the limit of 
small a. This merely reflects the fact that a truncated 
multipole expansion is valid only for slowly varying en- 
velopes. 

The corresponding results for Ino.53Gao.47As/InP are 
shown in Fig. 1201 Again, the error generated by the trun- 
cation of the multipole series vanishes away from the in- 
terfaces. The magnitude of the error is similar to that in 
GaAs/AlAs. For Ino.53Gao.47As/InP, however, there is a 
slight error in the bulk macroscopic electric fields even 
before the multipole expansion is truncated. 

To test the convergence of the truncated multipole ex- 
pansions used here, the power series (18b[) was extended 
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FIG. 19: (Color online) Errors in the macroscopic (a) electron 
density and (b) Hartree potential of a GaAs/AlAs superlattice 
under various approximations. The functions shown are the 
errors in the complete quadratic approximation to the macro- 
scopic density and potential (i.e., the macroscopic averages of 
the quadratic functions in Fig. [16] with a — 3d) and the corre- 
sponding errors when the linear and quadratic responses are 
generated from truncated multipole expansions. The labels 
< / < 4 and < I < 6 refer to the power series for the linear 
density; for the quadratic density, the corresponding limits 
are < I < 2 and < I < 4, respectively. 
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FIG. 20: (Color online) Errors in the macroscopic (a) electron 
density and (b) Hartree potential of an InP/Ino.53Gao.47As 
superlattice: comparison of complete quadratic errors with 
truncated multipole errors, as in Fig. 1191 



to include terms up to I = 6 in the linear density and up 
to I = 4 in the quadratic density. The resulting macro- 
scopic errors (shown in Figs.[T9land l2"0"|) actually increase 
slightly, but the change in the Hartree potential remains 
negligible. For the case a — 6d (not shown here), which 
includes the most important part of the potential for 
slowly varying envelopes, the macroscopic Hartree error 
of both truncated expansions is visually almost indistin- 
guishable from that of the complete quadratic response. 

As a direct practical test of the convergence of the 
multipole expansion, the extended limits (including also 
terms of the same order in the ionic pseudopotential and 
exchange-correlation potential) were used to recalculate 
the superlattice subband structures shown in Sec. V of 
the following paper^ No difference was observable in 
any case. Therefore, the truncated multipole expansions 
defined in Eq. (jlOp are very well converged for slowly 
varying envelopes. 



C. Bulk energy eigenvalues at F 

The error analysis to this point has focused on the 
Hartree potential, with the multipole truncation analy- 
sis limited to the macroscopic (G — 0) case. However, 
the Hamiltonian for slowly varying envelope functions is 
also sensitive to the G ^ terms, the quadratic-response 
error in the exchange-correlation potential, and the mul- 
tipole truncation error in this potential and the local and 
nonlocal ionic pseudopotentials. (There is no quadratic 
error in the latter because the ionic pseudopotential is 
purely linear.) 

These sources of error are not analyzed separately here, 
but to leading order their net effect is to generate an er- 
ror in the bulk energy bands of Figs. [2] and [3] This error 
can be calculated easily at the T point by adding linear 
and quadratic bulk perturbations to the T Hamiltonian 
of the reference crystal, diagonalizing the resulting matri- 
ces, and comparing the approximate T energies with the 
exact values for the given materials. For example, the 
quadratic error in the T^ v energy is —0.4 meV for GaAs 
and +0.5 meV for AlAs. Combining these values with 
the corresponding bulk macroscopic errors of +1.9 meV 
and —1.8 meV for the Hartree potential in Fig. [TWb). 
there is a net error of +1.5 meV in the bulk valence- 
band edge of GaAs and a net error of 2.7 meV (or 0.6%) 
in the GaAs/ AlAs valence-band offset. 

Likewise, the quadratic error in the Ti§ v energy is 
+3.15 meV for Ino.53Gao.47As and —4.45 meV for InP. 
The nominal positions of the interfaces in Fig. [501 are 
z = —0.25 and z = 23.75. The errors in the macro- 
scopic Hartree potential at the positions halfway be- 
tween the interfaces are +1.35 meV for Ino.53Gao.47As 
and —1.32 meV for InP. Therefore, there is a net 
quadratic error of +4.5 meV in the bulk valence-band 
edge of Ino.53Gao.47As and a net error of 10.3 meV 
(or 4.4%) in the Ino.53Gao.47As/InP valence-band offset. 
Note that the quadratic error in each component of the 
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Ino.53Gao.47As/lnP valence-band offset is only about 1%, 
but the total quadratic error is larger because the Hartree 
shift and bulk energy shift have opposite signs, while the 
corresponding errors have the same signs. 



VIII. CONCLUSIONS 

This paper has investigated numerically a quadratic- 
response approximation 10 ' 11 for the self-consistent one- 
electron potential within a model system based on su- 
perlattice LDA calculations with norm-conserving pseu- 
dopotentials. The electron density and potential energy 
of the superlattice were approximated by retaining only 
the linear and quadratic response to the heterostructure 
perturbation. This approximation worked well, with a 
net absolute error for the Tis valence states of about 2 
meV in GaAs/AlAs and 5 meV in Ino.53Gao.47As/InP. 
As shown in the following paper^ the principal effect of 
this error for slowly varying envelope functions is simply 
a constant shift of the superlattice energy eigenvalues. 

The density and short-range potentials were then ap- 
proximated further using truncated multipole expansions 
(i.e., power series in fc), retaining terms of order k 2 in the 
linear potential and k° in the quadratic potential. This 
had no effect on the macroscopic density and potential 
in bulk, but it generated some additional error (due pri- 
marily to the truncation of the linear density response) 
in a narrow region near the interfaces. This error was 
confirmed to be negligible for slowly varying envelope 
functions. 

Although the quadratic response is about an order of 
magnitude smaller than the linear response, it must be 
included if the self-consistent heterostructure potential is 
to be predicted to an accuracy of better than a few tens 
of meV. The quadratic response is also qualitatively im- 
portant because of its different symmetry. Dipole terms 
in the quadratic response were found to produce inter- 
face asymmetry and macroscopic electric fields in the no- 
common-atom Ino.53Gao.47As/InP system. As shown in 
the following paper ^ these terms, which have Ci v sym- 
metry, produce a significant fraction of the splitting of the 
quasidegenerate ground state in such systems. Therefore, 
to provide a realistic description of such interface-induced 
symmetry-breaking effects, existing empirical pseudopo- 
tentials would need to be modified to include these con- 
tributions. 



IX. LIMITATIONS AND POSSIBLE 
EXTENSIONS OF THE THEORY 

A. Quasiparticle effects 

The LDA model system used in the present study is, 
of course, not completely realistic because it fails to pro- 
vide accurate predictions of the conduction-band states. 



This model could be improved by including quasipar- 
ticle self-energies 3 -^ and by replacing norm-conserving 
pseudopotentials with projector augmented wave o 70 ' 71 ' 72 
(which would provide a better description of core elec- 
trons while remaining within a pseudopotential-like for- 
malism needed for the applicability of the present per- 
turbation scheme). Despite the well-known inaccura- 
cies of LDA conduction bands, 6 LDA wave functions 
do have a very high overlap with GW quasiparticle 
wave functions^ This similarity of the wave functions 
was used by Wang and Zunger— to construct empirical 
pseudopotentials with accurate energy gaps and effective 
masses using only a small (first-order) perturbation of 
the LDA Hamiltonian. This strongly suggests that, even 
though a numerical implementation of quasiparticle self- 
energies would be substantially more complicated than 
the present LDA model, the basic qualitative conclu- 
sion of the present study — namely, that the quadratic 
response provides an accurate approximation of the one- 
electron heterostructure potential — would remain valid 
in this theory also. 



B. Atomic relaxation 

Another useful extension of the present study would be 
to allow the superlattice ions to relax to their equilibrium 
positions. For example, it is well known that, although 
Ino.53Gao.47As and InP have the same bulk lattice con- 
stant, the bond lengths for the InAs or Ino.53Gao.47P 
bonds that form at a heterojunction are significantly 
different 46 ' 66 ' 74 The resulting displacement of the ionic 
planes from their present "clamped" reference-crystal po- 
sitions would generate additional dipole, quadrupole, and 
higher-order moments that have not been included here. 
The lowest-order contribution is a dipole term arising 
from the displacement of a single plane of ions, which 
generates a potential shif t 46 ' 75 AV = 4irZ*u/Ae 00 , where 
Z* is the transverse effective ionic charge, u is the dis- 
placement, A is the area per ion in the plane, and 
is the static electronic dielectric constant. The qualita- 
tive effects of such strain-induced interface dipoles are 
the same as those of the purely "chemical" dipoles in the 
clamped-ion system studied here, but accurate predic- 
tions of the properties of physical heterostructures could 
only be achieved by including both contributions. 

Atomic relaxation is readily calculated using any of 
several structural optimization algorithms already in- 
cluded in the ABINIT software package used here2&& 
The main difficulty in handling strained superlattices is 
the same as that encountered in strained bulk crystals — 
namely, that the electronic boundary conditions for the 
strained system are not the same as those for the un- 
strained system, so that the strain cannot be treated 
by the ordinary methods of perturbation theory i 76 ' 77 
This difficulty can be surmounted by an extensio n 78 ' 79 of 
the usual bulk coordinate-transformation technique 7 - 6 - in 
which the coordinates of one (or possibly more than one) 
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atom per unit cell in the strained system are mapped onto 
the coordinates of the same atom in the reference crys- 
tal. The coordinate transformation is straightforward in 
principle but carries with it an additional heavy layer 
of mathematical formalism. Therefore, in the interests 
of simplicity, this complication was excluded by fiat in 
the present work, although it may be taken up in future 
studies. 



C. Quantum wires or dots 

The final extension discussed here is the application 
of the present methods to quantum wires or quantum 
dots. This paper was limited to superlattices primarily 
for practical reasons, since, as mentioned in Sec. IIII Al 
these calculations are intended to be used in the follow- 
ing paper— to compare envelope-function predictions of 
valence subband structure directly with numerical LDA 
results. Such a direct comparison is possible for super- 
lattices containing on the order of 100 atoms, which are 
large enough for slowly varying envelopes to exist and 
yet small enough for LDA plane-wave calculations to be 
feasible. However, this clearly would not be practicable 
(at present) for a large quantum dot containing ~ 10 6 
atoms. 

Nevertheless, even if direct comparisons are currently 
out of reach for large quantum dots and wires, the 
quadratic-response approximation can still be used to 
construct envelope-function Hamiltonians for such sys- 
tems. Inspection of Fig. [TS] shows that the response to 
monatomic and diatomic perturbations is essentially in- 
dependent of the size of the supercell as long as the su- 
percell is large enough for interactions between perturba- 
tions in adjacent cells to be negligible. Therefore, one can 
construct the Hamiltonian from calculations on relatively 
small supercells containing one or two atomic perturba- 
tions (as opposed to the perturbations comprising one or 
two planes of atoms in the present study). Indeed, such 
calculations have already been performed for monatomic 
perturbations in 16-atom fee supercells in Refs. |45| and 
46, although somewhat larger supercells would probably 
be necessary for an accurate treatment of diatomic per- 
turbations. 

It should be noted, however, that the one-dimensional 
multipole expansions developed in Sec. I VII for superlat- 
tices are not directly applicable in quantum wires or dots. 
For such cases, one would need to use the more gen- 



eral three-dimensional multipole expansions developed in 
Refs. [l3 and [ill. These would lead to additional terms in 
the envelope-function Hamiltonian, as discussed in Ref. 
UTI and in footnote 30 of the following paper 
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APPENDIX 

If we consider for simplicity a function of only two 
variables, the density and potential are written as 
quadratic functions 

f{6 ll e 2 ) = c+c 1 e 1 +c 2 6 2 +c 11 e 2 l +c 22 el+2c 12 e 1 e 2 . (a.i) 

Here the coefficients c, c,-, and c^- were calculated by the 
direct supercell method^ First, the self-consistent den- 
sity and potential were found for a supercell consisting 
of the reference crystal with a monatomic perturbation 
6i = ±Aj (e.g., in Alo.5Gao.5As, one ALj.sGao.s atom was 
replaced by either Alo.55Gao.45 or Alo.45Gao.55). This 
provides the values of the coefficients 



ci 



Cll 



/oo, 
fio 



h 



2Ai ' 
/10 + /lo — 2/oo 
2A? 



(A.2a) 
(A.2b) 

(A.2c) 



where /j = /(— Ai,0), etc. Next, diatomic perturba- 
tions were used to find the value of c\ 2 . Two possible 
formulas were considered: 



C12 



C12 



./11 ~ hi - hi + hi 
8A!A 2 

fn — fio — foi + /oo 
2A1A2 



(A.3a) 
(A.3b) 



The symmetric formula (|A.3a[) is presumably more accu- 
rate, but very little difference was found in comparison to 
the asymmetric formula (|A.3b|) . Since the latter requires 
only one-fourth the computation time and storage of the 
former, Eq. (|A.3b|) was used for all calculations reported 
here. 
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